---
title: "Experiment 5a"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(tidyverse)
```

```{r import}
fulldata <- read_csv("~/Desktop/Studies/Force/analyses/FINAL/Exp.5/exp5a_completefinaldata_foranalysis.csv", 
    col_names = FALSE)


```
##Clean data
```{r clean_data}
#Name the columns
names(fulldata)<-c("sub","duration","captcha","scenario","normalt_abalt","alt1_norm","alt2_norm","alt1_val","alt2_val","force","comp1","comp2","sub2","scenario2","actual_act","alt1","alt2","comp3_resp","comp3")

#Remove rows with any missing data
cleandata<-fulldata[complete.cases(fulldata), ]

#Specify that the variables are factors
cleandata$scenario<-factor(cleandata$scenario,labels=c("heat","walks","weights"))
cleandata$sub<-factor(cleandata$sub)

#separate out scenarios for other analyses
heatdata=cleandata[which(cleandata$scenario=="heat"),]
walksdata=cleandata[which(cleandata$scenario=="walks"),]
weightsdata=cleandata[which(cleandata$scenario=="weights"),]

#create new alt variable that is the lowest of the two alt norm responses and highest of the two alt val responses
cleandata$critalt_norm<-ifelse(cleandata$alt1_norm < cleandata$alt2_norm,cleandata$alt1_norm,cleandata$alt2_norm)

cleandata$critalt_val<-ifelse(cleandata$alt1_val > cleandata$alt2_val,cleandata$alt1_val,cleandata$alt2_val)

#Create categorical norm and val data for future visualization
val_summary<-mosaic::favstats(~critalt_val,data=cleandata)
norm_summary<-mosaic::favstats(~critalt_norm,data=cleandata)

cleandata$val_cat<-ifelse(cleandata$critalt_val >= val_summary$Q3,cleandata$val_cat<-"high",ifelse(cleandata$critalt_val <= val_summary$Q1, cleandata$val_cat<-"low",cleandata$val_cat<-"medium"))

cleandata$norm_cat<-ifelse(cleandata$critalt_norm >= norm_summary$Q3,cleandata$norm_cat<-"high",ifelse(cleandata$critalt_norm <= norm_summary$Q1, cleandata$norm_cat<-"low",cleandata$norm_cat<-"medium"))

```
##Descriptives
```{r descriptives}

#get descriptives

#averaging across scenarios
force_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_avg)

norm1_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(critalt_norm),
               mean = mean(critalt_norm),
               sd   = sd(critalt_norm),
               se   = sd / sqrt(N)
)
print(norm1_desc_avg)

value1_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(critalt_val),
               mean = mean(critalt_val),
               sd   = sd(critalt_val),
               se   = sd / sqrt(N)
)
print(value1_desc_avg)

#between scenarios
force_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc)

norm1_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(critalt_norm),
               mean = mean(critalt_norm),
               sd   = sd(critalt_norm),
               se   = sd / sqrt(N)
)
print(norm1_desc)

value1_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(critalt_val),
               mean = mean(critalt_val),
               sd   = sd(critalt_val),
               se   = sd / sqrt(N)
)
print(value1_desc)

#with categorical data
force_desc_cat <- ddply(cleandata, c("val_cat","norm_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_cat)

force_desc_valcat <- ddply(cleandata, c("val_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_valcat)

force_desc_normcat <- ddply(cleandata, c("norm_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_normcat)

```

#Final analyses
```{r linear model of condition}

#Model predicting force from condition with random intercept for scenario
fullmod<-lme4::lmer(force ~ normalt_abalt + (1|scenario), data=cleandata)
summary(fullmod)

#Model excluding condition
modnocond<-lme4::lmer(force ~ (1|scenario), data=cleandata)
summary(modnocond)

#Compare the two
anova(fullmod, modnocond)

```

```{r linear model crit alternatives}

# Run full model with the critical alternative actions and random intercept for scenario
baseModel1 <- lme4::lmer(force ~ critalt_norm + critalt_val + (1|scenario), data=cleandata )
summary(baseModel1)

##Test the effect of value by dropping critalt_val from model and comparing to the original model:
modelNoVal1 <- lme4::lmer(force ~ critalt_norm + (1|scenario), data=cleandata )
anova(baseModel1,modelNoVal1)

##Test the effect of unusualness by dropping critalt_norm from model and comparing to the original model:
modelNoNorm1 <- lme4::lmer(force ~ critalt_val + (1|scenario), data=cleandata )
anova(baseModel1,modelNoNorm1)

#Add condition into the original model and compare to original 
baseModel_withcond <- lme4::lmer(force ~ critalt_norm + critalt_val + normalt_abalt + (1|scenario), data=cleandata )

anova(baseModel1,baseModel_withcond)
```
